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Abstract. A new numerical technique to identify the cosmic web is proposed. It is based 
on locating multi-stream flows, i.e. the places where the velocity field is multi- valued. The 
method is local in Eulerian space, simple and computaionally efficient. This technique uses 
the velocities of particles and thus takes into account the dynamical information. This is in 
contrast with the majority of standard methods that use the coordinates of particles only. 
Two quantities are computed in every mesh cell: the mean and variance of the velocity field. 
In the cells where the velocity is single-valued the variance must be equal to zero exactly, 
therefore the cells with non-zero variance are identified as multi-stream flows. The technique 
has been tested in a N-body simulation of the ACDM model. The preliminary analysis has 
shown that numerical noise does not pose a significant problem. The web identified by the 
new method has been compared with the web identified by the standard technique using 
only the particle coordinates. The comparison has shown overall similarity of two webs 
as expected, however they by no means are identical. For example, the isocontours of the 
corresponding fields have significantly different shapes and some density peaks of similar 
heights exhibit significant differences in the velocity variance and vice versa. This suggest 
that the density and velocity variance have a significant degree of independence. The shape of 
the two-dimensional pdf of density and velocity variance confirms this proposition. Thus, we 
conclude that the dynamical information probed by this technique introduces an additional 
dimension into analysis of the web. 
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1 Introduction 

The morphology of the cosmic web is an active field of research in modern cosmology. A wide 
variety of methods have been suggested, see e.g. [1, 2, 6, 7, 14, 17, 24, 26, 28] just to mention 
a few. However, all these methods share one common feature: they study the cosmic web 
as a property of the density field derived either from the distribution of particles in space in 
cosmological N-body simulations or from galaxies in the redshift space. However, the web 
can be viewed as a giant multi-stream flow. We suggest a simple technique that is capable 
to identify multi-stream flows and test it in a cosmological simulation of the ACDM model. 

We begin with a short overview of the origin and evolution of the concept of multi- 
stream flows in the context of the formation of the large-scale structure in the universe. The 
notion of three-stream flows was introduced by Zeldovich [32], in his paper on pancakes. He 
briefly discussed the formation of three-stream flows and illustrated the formation of a three- 
stream flow by the evolution of Eulerian coordinate as a function of Lagrangian coordinate 
with time. However, he himself used it mainly as a useful auxiliary theoretical tool helping 
to describe the nonlinear structures known at present as "Zeldovich's pancakes". Before 1980 
he himself seemed not to think that the multi-stream flows as a physical phenomenon might 
have relevance to the formation of the structures in the universe. Similar views were common 
in the west as well. For instance, multi-stream flows were barely mentioned in Peebles' book 
published in 1980 and it was only in the context of critical remarks on the pancake model. 

Doroshkevich et al. [12] explored the evolution of the multi-stream flows in the one- 
dimensional N-body simulation a little beyond the three-stream flow stage up to seven 
streams. This simulation has showen that the multi-stream flow region remains relatively thin 
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in the comoving coordinates in contrast with a simple extrapolation of the Zeldovich approx- 
imation (ZA) where it grows unlimitedly. The quantitative comparison of the thicknesses of 
the three-stream flow in ZA and in the one-dimensional numerical simulation was presented 
in the review [29]. The geometry of generic caustics in ZA in two dimensions was discussed 
in [3], the authors also provided a table of generic singularities occurring in the potential 
flows in three-dimensional space. A detailed study of the phase space in self-similar gravi- 
tational collapse in one, two and three dimensions was presented in [13]. A high resolution 
two-dimensional distribution of particles obtained in ZA that clearly showed the structure 
of multi-stream flows was demonstrated in [8], and in the high-resolution two-dimensional 
N-body simulations in [22]. 

The multi-stream flows were directly addressed in the model based on the adhesion 
approximation (AA) suggested in [15]. The AA model was designed to control the runaway 
growth of the thickness of the pancakes in ZA by introducing artificial viscosity term in the 
Euler equation. This term has no effect on the motion in voids leaving it as it was in ZA but 
does not allow the formation of multi-stream flows by annihilating the transverse to pancakes 
or filaments velocities. Thus, the high density walls and filaments are formed in AA instead of 
multi-stream flows in ZA. The density profile controlled by the adopted value of the viscosity 
coefficient becomes smooth and the velocity remains single valued. It is worth stressing 
that AA modifies the multi-stream flows by transforming them in single-stream flows of a 
special kind, characterized by high density. An important feature of the model consists in 
allowing the mass flow within the pancakes and filaments. The model incorporates most of 
the features of the hierarchical clustering process which is characteristic of the cosmological 
models dominated by cold dark matter. In addition, the AA model predicted several new 
features including the continuous flow of mass from walls to filaments to clumps, the multiple 
merger of clumps, the collapse of some voids, the presence of substructure of hierarchical 
type in voids, formation of the next generations of the filaments and pancakes. It was 
quantitatively tested against the two-dimensional [19] and three-dimensional [23] N-body 
simulations and was also used for predicting the structure in the forthcoming SDSS [30, 31]. 

Another significantly simpler modification of ZA was suggested in [9]. Since the ZA 
model has a serious problem with the runaway growth of the three-stream flow regions the 
authors proposed to filter out the perturbations with the scales smaller than the scale of 
nonlinearity i.e. the scale corresponding to the r.m.s density contrast fluctuation being 
equal to unity. The tests of the model dubbed the truncated Zeldovich approximation (TZA) 
against three-dimensional N-body simulations generally confirmed what was expected: the 
large-scale structure was reproduced quite accurately but the structures on small scales were 
erased. Two years earlier it had been shown that even replacing the small-scale perturbations 
with a new statistically independent realization would not change significantly the large- 
scale structure [20]. The large-scale structure proved to be quite robust. The evolution 
of the structure with time could be crudely probed by generating a sequence of particle 
distributions using TZA with the scale of nonlinearity corresponding to each chosen epoch. 

A very similar modification of ZA as far as the large-scale dynamics is concerned was 
suggested by Bond, Kofman and Pogosyan [5] (the BKP model). The only but crucial differ- 
ence between the BKP and TZA models consists in a different choice of the scale separating 
the large-scale and small-scale dynamics. The authors of BKP model strongly stressed that 
the boundary must be chosen in such a way that the large-scale motion was strictly single- 
stream flow. The authors also emphasized the coherence of the filtered linear density field 
and derived from it the strain tensor field on the scales of tens of Mpc between rare high 
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density peaks, the idea somewhat similar to the proposition made earlier in [10, 11] and [31]. 
The presence of this coherence causes the appearance of filaments and possibly pancakes. 
However, the chosen condition for the scale separating large-scale dynamics from small scale 
dynamics means that the filaments are the enhancements of density but not multi-stream 
flows. This is contrary to both the TZA and AA models. 

It is probably worth mentioning one more difference between AA on the one hand and 
TZA and BKP on the other although it is of a rather technical kind. The TZA and BKP 
models similarly to ZA use the the strain tensor as a primary initial field (referred to as 
the deformation tensor in TZA). The AA model is based on the linear velocity potential 
defined by the relation v = — V$ t) , which coincides up to a constant factor with the linear 
gravitational potential in the growing mode. In principle, both the velocity potential and 
strain fields contain exactly same information since one of them can be easily computed from 
the other. However, repackaging of the same information may be useful if it helps to identify 
the most significant variable that determines the structure or serves a particular goal of the 
model better. Summarizing this short discussion of three theoretical models of the large- 
scale structure we stress one difference essential for this paper. The multi-stream flows are 
characteristic for both TZA and AA while the BKP model eludes them in the large-scale 
dynamics. 

The ZA model defines pancakes as the regions between the surfaces on which the density 
is formally singular i.e. as the multi-stream flows. The three-dimensional N-body simulation 
by Klypin and Shandarin [18] revealed that the most conspicuous features after clumps are 
filaments rather than pancakes. Combining this finding with the results of [3] naturally 
caused the definition of filaments to become similar to that of pancakes: filaments are the 
multi-stream flows having very oblong shapes. The similarity is based on dynamics, the 
both pancakes and filaments are nonlinear but unvirialized concentrations of mass therefore 
in the dynamical hierarchy they hold an intermediate position between the density peaks 
that have not collapsed yet and virialized halos. At the same time another commonly used 
definition of objects as peaks above certain height in the density field [4] was also broadened 
to include pancakes and filaments [5]. On the one hand the former is more physical but 
at the time seemed to be harder to implement especially when the dynamical information 
about the structure was scarce. On the other the latter seems to be more practical but 
involves a free parameter, the threshold that determines which peak should be considered 
an object. This difference between two definitions was probably the reason of the strongest 
discrepancy between the models declared in [5]: the formation of the objects in BKS is 
in inverse order than in ZA. It is worth pointing out that despite considerable overlapping 
between two definitions the objects they describe are not identical. 

Despite the inability of current N-body simulations to resolve properly the phase space 
there is no doubt that the multi-stream flows must be present in the collisionless medium 
in the nonlinear regime. Therefore even the limited information about multi-stream flows 
would introduce a new dimension into the analysis of the cosmic web as it does in virialized 
halos [21]. This must be based on using the velocity field. Velocities bring about dynamical 
information totally independent of the density and gravitational potential fields in a general 
case. There are of course special cases like the systems in virial or thermal equilibrium 
where certain relations of statistical nature between velocities, coordinates of the particles 
and gravitational potential could be derived. 

Recently there have already been made some attempts to incorporate dynamical infor- 
mation into the analysis of cosmic web [14, 17]. Both approaches are based on the analysis 
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of the eigen values of the Hessian of the gravitational potential at the nonlinear stage. Both 
groups appeal to the analogy with ZA and claim that this statistics provides a dynamical 
classification of the cosmic web. The only difference between two approaches is in the dif- 
ferent choices of the amplitude to be assigned to the eigen values of the Hessian. However, 
as we have already mentioned neither density nor gravitational potential contain dynamic 
information in a general case since the both are computed from the particle coordinates only 
and absolutely independent of the velocities. We have already mentioned a couple of special 
cases when the gravitational potential and velocities could be related in statistical sense. The 
growing mode in the linear regime represents another special case. The dynamical evolution 
in ZA is described by a map x(q, t) = q — D(t)V^ v (ci), where $> v is the velocity poten- 
tial not the gravitational potential. It is only due to the degeneracy of the growing mode 
in the linear regime where the velocity potential proportional to the gravitational potential 
& v = constant x & g they can be used interchangeably with the proper choice of the constant. 
Therefore the gradient of the gravitational potential is proportional to the velocity and it 
can be computed from the gravitational potential in the growing mode in the linear regime. 
The ZA model is of course an extrapolation and therefore required additional scrutiny before 
it was excepted as a reasonable approximation for the initial phase of the nonlinear regime 
(see e.g. [29]). However, neither the gradient of the potential can be used for mapping nor 
its second derivatives represent the deformation tensor in a general case. If correct the rea- 
soning presented in [17] and [14] must be valid for arbitrary velocities of the particles since 
the analysis and thus the conclusions do not depend on the velocities at all. This hardly 
can be true. Thus, even the suggested statistics may provide some merits for the study of 
the density and gravitational potential fields at the nonlinear stages its interpretation as a 
dynamical characteristics cannot be excepted in the present form. 

The goal of the paper is to introduce a new numerical technique to identify the cosmic 
web as the multi-stream flows and present the results of the first tests, which are very 
encouraging. The method is based on a simple statistics derived from the coordinates and 
velocities of the particles: the mean and variance of velocities in every mesh cell. Here we 
use a uniform spatial mesh although the homogeneity of the mesh is not required by the 
method. We estimate the role of noise caused by numerical effects and compute the two- 
dimensional pdf of the density and variance of velocity which demonstrates a significant level 
of independence of these quantities. We compare the appearance of the multi-stream flow 
web with the web derived from the density field only as well as with the distribution of the 
particles themselves. 

Although there is no agreement between different groups on the exact definition of the 
filaments and walls/panckes the most agree that the filaments comprise more mass than any 
of the rest constituents of the web (i.e. clusters /clumps, walls/pancakes, field/void). For 
instance, the fraction of mass in filaments and walls in the ACDM simulations to be almost 
50% of the total mass, considerably more than in clusters [1]. Both filaments and walls are 
nonlinear but unvirialized objects in contrast to clusters/clumps which are mostly virialized 
objects. The suggested method allows to quantitatively analyze the dynamics of these objects 
because they are multi-stream flows. This work is under way and the results will be reported 
in detail in a separate paper. 

The rest of the paper is organized as follows. Section 2 describes the idea of the method, 
Sec. 3 describes the numerical technique and N-body simulation used in the tests, the noise 
analysis is given in Sec. 4. The total volume devoid of multi-stream flows is evaluated in 
Sec. 5. The appearance of the multi-flow web is compared with the density based web in 
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Figure 1. Top left panel. The density as a function of the Eulerian coordinate, x, is shown for all 
streams by dotted and dashed lines. The total density is shown by a solid line. Top right panel. 
The velocity in three streams are shown by dashed and dotted lines, the solid line shows the velocity 
outside the three-stream flow. Bottom left panel. The solid and dashed lines show respectively 
and e defined in eq. 2.2. Bottom right panel shows the mean velocity defined in eq. 2.3. 

Sec. 6, and the mean velocity field is shown in Sec. 7. Finally the results are summarized in 
Sec. 8. In the rest of the paper density is given in the units of the mean dark matter density, 
velocities in km/s and a 2 , in (km/s) 2 . 

2 The idea of the method 

The physical idea of the method can be illustrated by a simple example. Let us consider the 
formation of a one-dimensional pancake in ZA for a simple sinusoidal perturbation [3, 29, 
32]. Using the comoving coordinates one can easily write down explicit expressions for the 
position, velocity and density as a function of the Lagrangian coordinate q and the linear 
growth factor D that can effectively play the role of time 

x = q — Dsinq, v = dx/dD = — sing, p = |1 — D cos q\~ ■ (2-1) 

The solution predicts a shell-crossing singularity in density at x = at D = 1 and the 
formation of a three-stream flow afterward. At small 6D = D — 1 it approximately describes 
the generic phase space structure of the shell crossing regions at their early stages. At later 
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times it evolves into a multi-stream flow with five then seven etc. streams which is beyond 
the scope of ZA. The density and velocity profiles in the three-stream flow and its vicinity 
are shown in the top panels in Fig. 1. The velocity and density are shown for each stream in 
addition the total density in the three-stream flow region is also shown. Two bottom panels 
show the mean velocity, velocity variance a 1 (solid line) and the density of kinetic energy e 
(dashed line) defined bellow. 

Usually the analysis of the structure of multi-stream flows is limited to the discussion 
of velocity and density fields. We propose to examine two additional quantities of thermo- 
dynamic and therefore statistical kind similar to density. One of them is the variance of 
velocity, cr^(x), and the other is the mean kinetic energy e(x) at every Eulerian point x 

where Pi(x) and i>i(x) are the density and velocity in each stream and Avi = V{(x) — v(x), 
n s is the number of streams. The mean velocity is defined as 

Hx) s S&ffM^M . (2.3) 

Z_/i=l Pi\ x ) 

Figure 1 shows all these quantities at the nonlinear stage when the shell crossing has already 
happened and the three-stream flow has formed. Passing by we note that the mean velocity 
of the three-stream flow does not coincide with the velocity in one of the streams as may 
appear in the figure, but they are quite similar in this particular example. 

Comparing the shapes of <r^ and e curves with the density and velocity curves we see 
remarkable differences. The most important, however not surprising feature of the and 
e curves is that they are identical zeros beyond the multi-stream flow regions. Less obvious 
and more surprising feature is that their shapes are so similar to each other and practically 
inverse to the shape of the total density curve in the shell crossing region. The both features 
suggest that the and e functions could be complimentary characteristics to the density 
in the studies of the dark matter cosmic web. The fact that they are equal exactly to zero 
beyond the multi-stream flow regions may mean that they can be good tracers of multi- 
stream flows. However, the discussed example is too oversimplified. In the reality the noise 
caused by discreteness undoubtedly violates this ideal situation and therefore must be taken 
into account before one can come to a sound conclusion. But first we describe the numerical 
technique. 



3 The numerical technique and N-body simulation 
3.1 Numerical technique 

Since a% and e are so similar we shall limit the discussion to a% only. The velocity variance 
depends on density less than e and therefore it carries more independent information. The 
physical dimensions of a% suggest that it can be directly related to the gravitational potential. 
However a more thorough analysis may reveal some additional attractive features of e. 

We use the standard CIC (cloud-in-cell) method for the evaluation of v and then <r^. 
The particles are modeled as constant density cubes of size I = L/N p where L =512 7i is 
the size of the simulation box and N p =512 is the number of particles in one dimension. 
Both v and a?, are evaluated on a uniform cubic mesh with the cells of the same size as 
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that of particles. However, it is worth mentioning that the uniformity of the mesh is not a 
requirement of the method. In general each particle contributes a volume weighted fraction 
of its velocity to eight neighboring mesh sites. Each velocity fraction equals the fraction of 
the volume of the overlap of the particle cloud with the mesh cell. It also can be viewed as 
a fraction of its linear momentum since all the particles have the same masses. The mean 
velocity assigned to the mesh point is the sum of the contributions from all the particles 
overlapping with it divided by the total mass in the cell. The variance a 2 is evaluated in a 
similar manner. It is worth stressing that we conduct the further analysis without additional 
filtering of the density field which is also computed using the CIC scheme. 

3.2 The cosmological model and simulation parameters 

We applied this method to the pure dark matter N-body simulation in the 512 h M.pc 
cubic box using PM (particle mesh) code [16, 27]. The number of particles was 512 3 and the 
mesh in the gravitational force solver was 1024 3 . The parameters of the ACDM cosmological 
model were as follows: h = F /(100 km/s ■ Mpc) = 0.72, fl tot = 0.25, n b = 0.043, n = 0.97, 
<7g = 0.8, the initial redshift z- m = 200. The choice of the parameters is related to the main 
purpose of the study: to test this new technique and illustrate its performance by applying 
it to a realistic cosmological model. 

4 The analysis of numerical noise 

We begin with applying this method to the initial state when all the fields are in the linear 
regime. 

4.1 Initial linear state 

Even at the very beginning the small displacements of particles from unperturbed positions 
on the regular grid result in small overlapping of particles with each other. It is also easy to 
imagine that more than one particle may overlap a mesh cell. This results in generating of 
small but nonetheless nonzero a 2 . The left hand side panel in Fig. 2 shows three functions 
f(Vp) = n(t>p)/512 3 (in blue), f(v 2 ) (in green) and f(a 2 ) (in red), which are the fractions of 
particles or cells per equally sized logarithmic bins. In the ideal situation one expects that 
f(vp) = f(v 2 ) and f(a 2 ) = as it was described in Sec. 2. The first condition is satisfied 
with high accuracy and the second is obviously violated. However, the a 2 is approximately 
three orders of magnitude smaller than v 2 or v 2 therefore we conclude that this test is passed 
without serious problem. In the multi-stream flows the a 2 is expected to be less than v 2 but 
not orders of magnitude smaller (see Fig. 1 for an illustration). 

An additional question may occur. Does the generation of noise in a 2 depend on the 
density in the cell? Naively one may think that a 2 is generated in the regions where particles 
are crowded and much less in the underdense regions with p < 1 where the particles are 
deserted. The answer to this question happenes to be not that simple as the left hand 
side panel in Fig. 3 shows. The figure displays the two-dimensional histogram f(p,a 2 ) for 
the initial state. The maximum of the function is at p ~ 1, log 10 (<7^) ~ 5.5. The peak 
has approximately a triangle shape, extending to somewhat higher values of a 2 for both 
higher and lower densities. The third direction is downward to lower values of log 10 (cr^). 
This relatively complicated shape of f(p, a 2 ) indicates that the kinematics of particles is not 
simply contraction for p > 1 and expansion for p < 1 even at this early stage. The overall 
shape of the histogram shows no strong correlation between a 2 and p. 
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Figure 2. The fraction of particles per bin as a function of v 2 (blue histograms), the fraction of 
mesh cells per bin as a function of v (green histograms), and as a function of a 2 (red histograms) 
are shown for the initial state in the left hand side panel and for the present time in the right hand 
side panel. The fraction of particles with v 2 and the fraction of cells with v 2 since v 2 = v 2 in the 
one-stream flow regions, thus the green and blue curves are one on top of the other in the left panel. 
There are also 33.4% of cells with a 2 < 10 -9 not shown in the right hand panel. 
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Figure 3. The fraction of cells n(p, cr^)/512 3 in bins as a function of p and a 2 . The color bars 
show the logarithmic scales for each panel. The both histograms are of same size 50x50 with equally 
spaced bins along the axis marked in each panel. Left: the initial linear state. Right: the final highly 
nonlinear state. 

4.2 Final nonlinear state 

The statistics f{v 2 ), f{v 2 ), f{<T 2 ), and f(p,a 2 ) look quite differently at the nonlinear stage. 
The right hand side panel in Fig. 2 shows the histograms for the first three functions. 
Only f(Vp) and f(v 2 ) bellow the maximum are similar. At the high value end f(v 2 ) drops 
considerably faster than f(v 2 ) since the mean velocity in the multi-stream flow regions is 
significantly lower then the velocities of particles (see also Fig. 1). The high value end of 
f{<J 2 ) is much closer to that of both f(v 2 ) and f(v 2 ) as one may naturally expect for the 
multi-stream flow regions. It is worth mentioning that the highest mean velocities are greater 
than the highest relative velocities inside multi-stream flows. 

As many as 11.76% of cells have a 2 = exactly in excellent agreement with the fraction 
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Threshold 



Figure 4. The fraction of volume where er 2 , (solid line) or p/p (dashed line) is less than the threshold. 
The horizontal gives the threshold in (km/s) 2 for cr 2 and in dimensionless units for p/p. The low limit 
for both quantities is the same 11.7% which is the fraction of sites where the density is zero exactly 
in the N-body simulation. 

of completely empty cells with p = being equal to 11.74% , in addition there are 21.6% of 
cells with < a% < 1 x 10~ 9 . This explains the presence of red bar at the bottom of Fig. 3 
(right panel). The two-dimansional histogram f(p, <t^) shown in the right hand side panel of 
Fig. 3 is also substantially different from the linear case. There is an easilly seen correlation 
between and p, at p > 1. However, the spread is about two orders of magnitude along 
the both directions for the red peak. The former is naturally related to the correlation of 
the gravitational potential depth with the high density peaks. The latter suggests that many 
medium density peaks (10 < p < (a few) x 100) as well as the filaments are not relaxed 
systems with in such cells to be almost independent of p. 

5 The total volume in voids 

There are 11.74% of the mesh cells completely empty, they obviously have a% = exactly 
as well (except a tiny fraction due to numerical errors). However, additional 21.6% of cells 
have such a tiny a% < 10 -9 that they can be considered belonging to one-stream flow regions 
without doubts. The fractions of cells one as a function of density and the other as a function 
of where the corresponding quantity is less then the threshold shown in the horizontal 
is shown in Fig. 4. The function of a% grows from 33.24% at a% = 1 X 10~ 9 to 33.67% at 
a\, = 1 and to 35.26% at <r^ = 10. Then the rate of growth becomes significantly greater and 
the fraction f(a^) reaches 46.27% and 60.70% at =100 and 316 respectively. The long 
plateau from a% = 10~ 9 to about 10 (ten orders of magnitude!) shows that in about one 
third of the volume the velocity is a single valued function, i.e. no multi-stream flow regions 
are present there. This number can serve as a low physical limit for the total volume in 
voids for the resolution of this simulation. For a given mass resolution of N-body simulations 
one can give a physical definition of voids as the regions where the shell-crossing has not 
yet occurred. Therefore neither halos nor filaments nor sheets could have formed in these 
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Figure 5. A part of a thin slice of 1 h^ 1 Mpc thickness through N-body simulation box. The hllcd 
contours show the values of tr^ as indicated by the color column. Heavy red contours show the CIC 
density held with p = 1, the white thin contours inside red contours show p — 3, 10 and 50 The labels 
show the distances in h^ 1 Mpc. 

regions. The multi-stream flow regions can apparently appear in such places with the growth 
of the resolution of the simulation. 

6 The multi-stream flows versus the density web 

In this section we compare the multi-flow web with the standard representation of the web by 
the isocontours of the dark matter density. Figure 5 shows the density contours superimposed 
on the field of a\ in a thin slice (1 /i -1 Mpc) randomly selected from the simulation cube. 
The filled contours show the regions of the <r^ field between the levels marked on the color 
column. The density contours correspond to the following levels: red to p = 1 and three white 
contours correspond to p = 3, 10 and 50. One can see that although the overall structures are 
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very similar, but they by no means are identical. We consider the both observations as good 
news. The former is the evidence of the robustness of the identification of the filaments and 
other structures while the latter suggests that dynamical information in the a% field brings 
an additional dimension to the characteristics of the web. A closer inspection show that the 
highest peaks in two fields do not have identical shapes. It will be better seen in the following 
figures where we zoom up two regions: one with the lowest mean density (void region) and 
the other with the highest density in the slice (clumps and filaments). 

However, we first show the superposition of particles contributing to the density and 
in this slice on the a% field in Fig. 6. The color map is slightly different from the previous 
figure, which allows to see the structures at low values of better. In general, we conclude 
that the correspondence of two fields it very good as it should be apart from numerical noise. 
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Figure 7. A zoomed low density region in the top right corner of Fig. 5 and 6. The filled contours 
show the values of a% as indicated by the color column. The black and white contour lines show 
a few density peaks above unity in units of the mean density. The solid, dashed and dotted lines 
correspond to p = 0.5, 0.25 and 0.1 respectively. The heavy white line at about (130, 105) shows one 
small contour with p = 2. 

6.1 The slice through a void 

The region from the top right corner of the slice is zoomed in Fig. 7. The color map is 
changed again in order to adopt better for the low density environment. The levels of are 
shown by the color column while the density contours are shown by black lines. The density 
in the most of the region is below the mean except a few peaks shown by black and white 
heavy line. The impression is similar to that from Fig. 5: the fields have many similarities 
but also substantial differences, perhaps even a little greater differences than in Fig. 5. In 
particular, the empty regions seem to be significantly more empty in cr^ than in p in an 
agreement with the discussion in Sec. 5. Although the density contours are pushed to such a 
low values that one may completely disregard them as being subgrid noise the values of the 
<jy contours are also very low and still there is some similarity between both fields in some 
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Figure 8. A zoomed high density region in the bottom left corner of Fig. 5 and 6. Here, filled 
contours correspond to as indicated by the color. The black and white line shows contours of p 
= 1, and white solid, dashed, dotted and heavy solid contours correspond to p =3 , 10 , 30 , and 60 
respectively. 

parts of the figure. Thus, an optimistic conclusion would be that the combination of both 
fields may allow to probe low density regions more reliably than any of them alone. On the 
other hand, a pessimistic conclusion would relate it to the correlation of the p and fields 
at small a^. The further study is obviously required. 

6.2 The slice through clumps and filaments 

Figure 8 shows the relatively high density region in the bottom left corner of Fig. 5 and 6. The 
figure contains a few density peaks with p > 10 (white dashed lines) three of which are higher 
than p = 30 and two even higher with p > 60. There are also a few filaments connecting 
the clumps. Again the both fields are rather similar but there are notable differences as 
well. Two highest density peaks in the bottom left corner reach p > 60 but only one exceeds 
o"^ > 10 5 . However, the peak of extends to the density level as low as the mean density. 
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On the other hand a little lower density peak with p < 60 coincides with a quite substantial 
region with cr^ > 10 5 . Similar observations can be made with the peaks of lower density. 
Thus, we conclude that the density and are two parameters that have a substantial degree 
of independence and therefore may better characterize the clumps and filaments as any of 
them alone. In addition, the figure shows that the filaments are the multi-stream flows in 
agreement with both the TZA and AA models and contrary to the prediction of the BKP 
model that implies that the filaments have not reached the state of multi-stream flows. 

7 The mean velocity field 

Finally, we present the mean velocity field defined by eq. 2.3. We superimpose the mean 
velocity field shown by the arrows on Fig. 8 where we keep only three highest levels of the 
density contours for clarity. 

The highest density peak sitting at the largest patch of high is in the bottom left 
corner approximately at (17, 39). The overall pattern of the mean velocity field clearly 
indicates that the whole region shown in the figure is falling on this clump. The mass inside 
three filaments is streaming toward the clump confirming one of the predictions of the AA 
model [15]. The second largest clump at (28, 54) is falling on the first clump accompanied by 
the surrounding filaments. The filaments are highly inhomogeneous and the streams are far 
from being uniform. The velocity field in the vicinity of the first clump looks roughly circular 
giving the impression of quasi-spherical (actually quasi-circular in this plane because we see 
only the cross section) accretion but this is misleading if we look at the flux of mass pv shown 
by arrows in Fig. 10. It is worth stressing that in order to accommodate the arrow lengths in 
the figure, their lengths are made proportional to the square root of the magnitude instead 
of the magnitude. Thus the difference between the fluxes from three adjacent filaments 
compared to that from three surrounding voids is considerably more dramatic than it appears 
in the figure. 

8 Summary 

We present a new numerical technique that allows to identify the multi-stream flows, i.e. the 
regions where the velocity field is multi- valued. According to [32] these regions have reached a 
clear physical threshold to be identified as being in strong nonlinear regime but most of them 
are not virialized. The multi-stream flows form the web overall very similar by appearance 
to the web identified in the density field. At the same time many features are very different 
both in low and high density environments as seen in Fig. 7 and 8. These include the 
difference of the shapes of the contours around peaks, the lack of monotonic relation between 
the heights of peaks in the p and <t^ fields. Although one can see a positive correlation 
between the density and the a% that determines the multi-stream flows in Fig. 3, the spread 
of points along both axes (p and ) is very large suggesting that these parameters are quite 
independent. Therefore may serve as a second parameter characterizing the dynamical 
environment in the nonlinear density peaks before they reached virial equilibrium as well as 
for filaments and pancakes. The method for computing a% is localized in Eulerian coordinates 
and is very easy to implement. 

Ideally is identically equal to zero everywhere where the velocity is singled- valued but 
the simplest numerical implementation based on the ordinary CIC method generates noise. 
Although it is seemingly not posing a serious problem other schemes like TSC (triangular 
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Figure 9. The mean velocity field v projected on the slice plane is superimposed on the a% field 
shown in Fig. 8. White contours show peaks with density above 10, 30 and 60 in the units of the 
mean density with dashed, dotted and heavy solid lines respectively. 

shape density cloud) may be less noisy. This will be checked an reported in the following 
publications. 

The field naturally defines the physical condition for voids as the regions devoid of 
multi-stream flows. The identification of the multi-stream flows of course depends on the 
resolution (in particular the mass resolution) of the simulation, but this is the fundamental 
limitation for all fields obtained in the simulations. The a\, field by no means is worse than 
commonly used density field in this respect. This condition is physically more clear than one 
based on an essentially arbitrary density threshold. It is worth stressing that about 22% of the 
total volume has nonzero density but no multi-stream flows. For the resolution of the ACDM 
simulation used here (~1 /i~ 1 Mpc) we find the volume fraction in voids is at least about 
one third of the total volume. This is considerably greater than at the percolation transition 
evaluated for a similar simulation in [28], therefore the volume devoid of multi-stream flows 
consists predominantly of one connected region. 
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Figure 10. The flux of mass /?v is shown by arrows. The color scheme and contours are the same as 
in Fig. 9. In order to accommodate the length of arrows their lengths are made proportional to the 
square root of the magnitude. 

The mean velocity field Vi defined by eq. 2.3 shows the mass flow along the filaments 
toward the clumps in agreement with the prediction of the adhesion approximation [15]. 
Some filaments seem to move predominantly in the transverse direction as a whole (see Fig. 
9). The mean velocity field seems to be quite smooth without significance convergence in 
these cases. A somewhat similar observation can be made concerning the clumps. In some of 
them as in one at (17, 39) in Fig. 9 and 10 the velocity field seems to converge significantly 
stronger than in the other similar clump at (28, 54) in the same figures where the velocity 
field looks significantly more steady. 

The filaments with lengths over 50 /i _1 are shown to be the multi-stream flows in a 
qualitative agreement with both the truncated Zeldovich approximation [9] and the adhesion 
approximation [15]. 

Summarizing, we would like to stress the importance of the study of dynamics in fil- 
aments and pancakes, which are highly nonlinear but unvirialized dynamical systems. The 
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suggested technique may be a useful tool to serve this purpose. 
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